Problem Set 5: Model Selection

PLS 206 — Applied Multivariate Modeling

Author

Grey Monroe

Published

January 1, 2026

Note

Due: Monday, November 2, 1:10 PM (start of class) — submitted through Canvas.

Questions 1–4 use Monday’s lecture (AIC by hand, Akaike weights, AICc and BIC, testing a small set of hypotheses with model.sel()). Questions 5–10 use Wednesday’s (collinearity and VIF, dredge(), importance and model averaging, stepwise selection and its traps). The week’s code page has working code for every step; the hints point to its parts and sections.

Submission Instructions

Upload two files to Canvas:

Written answers — a .pdf document with your written responses, interpretations, and figures.
R script — a .R file with all the code you used to produce your answers.

List the first and last names of any collaborators in both files.

Naming convention: HW5emailID.R (code) and HW5emailID.pdf (written answers), where emailID is the part of your UC Davis email before the @.

Example: Grey Monroe would submit HW5gmonroe.R and HW5gmonroe.pdf.

Feel free to use Google, StackOverflow, and other resources if you are stuck with code, and post questions on Slack if anything is unclear.

Every model you plot: show the raw data, the fitted curve, and its 95% confidence band, with the response on the y-axis.

The data: flowering time in wild wheat (again)

You met these data in Problem Set 2. Flowering is flowering time (days from germination to flowering) for 54 genotypes of wild einkorn wheat (Triticum monococcum) collected across the Middle East and Europe, measured in a common garden: seeds from every collection site were grown together in a greenhouse, so differences among genotypes are genetic, not environmental. For each genotype’s collection site we have Latitude, Longitude and the 19 bioclimatic variables of its home environment from WorldClim (bio1 to bio19).

In Problem Set 2 you tested one climate variable at a time. This week the question is which explanation is best supported, and what happens when we let the computer choose.

library(ggplot2)
library(MuMIn)   # model.sel(), dredge(), sw(), model.avg()
library(car)     # vif()
library(MASS)    # stepAIC()
wheat <- read.csv("https://greymonroe.github.io/PLS_206/data/wild_wheat.csv")

Part 1: Three Hypotheses

Three ideas for what sets flowering time in wild wheat, each with one variable to stand for it:

Temperature: warmer sites select for faster development (bio1, annual mean temperature).
Summer drought: plants from sites with dry summers flower early to escape the drought (bio18, precipitation of the warmest quarter).
Day length: photoperiod changes with latitude, and flowering is often cued by day length (Latitude).

Your candidate set is six models: the null (Flowering ~ 1), each idea alone, temperature + drought (bio1 + bio18) and temperature + day length (bio1 + Latitude).

Q1. Fit the temperature model, lm(Flowering ~ bio1, data = wheat). Using logLik(), calculate its AIC by hand with AIC = −2 log L + 2k, and check your answer against AIC(). What is k for this model, and what does it count?

Hint

as.numeric(logLik(m)) gives log L, and attr(logLik(m), "df") gives the k that R uses. Count the parameters the model estimates: there is one more than the coefficients in coef(m). Code page Part B, sections 6 and 7; slides logLik() in R and AIC() does the same sum.

Q2. Fit all six models and compare them with model.sel(). Which hypothesis is best supported, and what is its Akaike weight? The two-predictor models are both about 2 AICc units behind the temperature model: compare their log-likelihoods with the temperature model’s and explain where the gap of about 2 comes from. Do drought or day length add anything once temperature is in the model?

Hint

model.sel(m0, m_temp, m_dry, m_day, m_temp_dry, m_temp_day). The logLik column is the fit, the df column is k. Slides Only differences matter: ΔAIC, Akaike weights: share of the support, model.sel(): the six models, ranked and The evidence, as weights; code page Part B, section 10, and Part C, sections 13 and 14.

Q3. Rank the same six models by AIC and by BIC as well as by AICc (model.sel(..., rank = AIC), rank = BIC). Does the ranking change? How does the weight of the temperature model change from AIC to AICc to BIC, and why? With n = 54, which of AIC and AICc should you report?

Hint

AICc = AIC + 2k(k + 1) / (n − k − 1); BIC = −2 log L + k log(n). Compare each one’s penalty for one extra parameter: log(54). Slides AICc: a correction for small samples, BIC: a stiffer price and AIC or BIC?; code page Part B, sections 11 and 12, and Part C, section 13 (rank = BIC).

Q4. Three comparisons. For each, say whether you can use anova(), AIC, both, or neither, and give the result where you can:

1. temperature alone vs temperature + drought;
1. temperature alone vs drought alone;
1. a labmate fits lm(log(Flowering) ~ bio1, data = wheat), finds an AIC of about −146 against 353 for your temperature model, and says the log model is far better.

Hint

anova(m_temp, m_temp_dry). anova() needs one model to be a special case of the other. For c), ask what the likelihood is a likelihood of. Slides Nested models: anova() and AIC agree, Not nested? Only AIC can compare them and Same response: log(y) and y cannot be compared; code page Part C, sections 15, 16 and 18.

Part 2: Exploring All 19 Climate Variables

Q5. Fit the model with all 19 climate variables. bio7 gets NA (you found why in Problem Set 2: bio7 = bio5 − bio6), so drop it. Then prune the other 18 with vif(): drop the variable with the highest VIF, refit, and repeat until every VIF is below 5. Which variables are left? Was bio1, the best single predictor of flowering time, kept? Why did the pruning throw it out?

Hint

A repeat loop saves you 12 rounds of copy-paste:

vars <- setdiff(paste0("bio", 1:19), "bio7")
repeat {
  v <- vif(lm(reformulate(vars, "Flowering"), data = wheat))
  if (max(v) < 5) break
  vars <- setdiff(vars, names(which.max(v)))
}

For the last part, look at cor(wheat$bio1, wheat[, c("bio5", "bio6", "bio10", "bio11")]). Slide Drop the worst VIF, refit, repeat; code page Part D, section 25.

Q6. Run dredge() on the model with the variables left after Q5. What is the top model by AICc, and which models are within 2 AICc units of it? Rerun with rank = "BIC": is the top model the same, and how does the runner-up differ? Plot Flowering against the top model’s predictor with the smallest p-value, with the top model’s fitted line and 95% band (hold the other predictors at their means).

Hint

options(na.action = "na.fail") first, or dredge() refuses to run. d <- dredge(global); subset(d, delta < 2); dredge(global, rank = "BIC"); top <- get.models(d, 1)[[1]]. For the band, build newdata with that predictor on a grid and the others at mean(), then predict(top, newdata, interval = "confidence"). Slides dredge(): every subset of the 7, The top 10 of 128 and AICc vs BIC with 6,299 records; code page Part D, sections 27 to 29, and Part B, section 6 for the band.

Q7. Report each variable’s importance from sw() on the full dredge table, and average the models within 2 AICc units (model.avg(d, subset = delta < 2)). Every variable appears in exactly half of the 64 models; what importance would a variable that explains nothing get, and what does that say about the importance of bio12? In the model average, why is the bio12 slope smaller in the full average than in the conditional average?

Hint

sw(d); summary(model.avg(d, subset = delta < 2)) prints both averages. The full average counts a slope of zero for every model that leaves the variable out. Slides Importance: summed weights and Model averaging: estimates across the top models; code page Part D, sections 30 and 31.

Q8. Your Part 1 temperature model and the top dredge model from Q6 were fit to the same 54 genotypes with the same response, so their AICc values can be compared. Which is better supported, and by how much? Is the top dredge model a finding or a hypothesis, and which of the three goals (exploration, inference, prediction) does the dredge in Q6 serve? Two or three sentences.

Hint

AICc(m_temp) and AICc(top) (both from MuMIn). Think about which variables dredge() was allowed to try. Slides A tool for each goal and What we found: hypotheses, not results; code page Part A, section 5, and Part D, section 32.

Part 3: The Traps

Q9. Stepwise selection follows one path. (a) Run stepAIC() backward from the 18-variable model (all bio variables except bio7), and forward from the null model with the 18-variable model as the upper scope. Do they agree? Compare their AIC values. (b) Now take the six variables from Q5, resample the 54 genotypes with replacement 200 times (set.seed(206) first), rerun backward stepAIC() on each resample, and tabulate how often each variable is kept. Which variables are stable, and how often does the most common model come up?

Hint

stepAIC(full18, trace = 0) goes backward; stepAIC(lm(Flowering ~ 1, data = wheat), scope = list(lower = ~ 1, upper = formula(full18)), direction = "forward", trace = 0) goes forward. Compare with AIC(), not the AIC column stepAIC() prints (it drops a constant). (b) In a loop: wb <- wheat[sample(nrow(wheat), replace = TRUE), ], refit and run stepAIC(), then attr(terms(s), "term.labels") gives the kept variables. Slides Stepwise: one path through the maze and Resample the farms, rerun stepAIC; code page Part E, sections 33 and 34.

Q10. Selection finds patterns in noise. Add 10 columns of pure noise to the six Q5 variables and draw a random half of the genotypes, running this block exactly as written:

set.seed(19)
wn <- wheat[, c("Flowering", "bio3", "bio4", "bio6", "bio8", "bio12", "bio18")]
for (j in 1:10) wn[[paste0("noise", j)]] <- rnorm(nrow(wn))
half <- sample(nrow(wn), 27)

Run backward stepAIC() on lm(Flowering ~ ., data = wn) using all 54 rows. Is any noise column kept with p < 0.01? (b) Now select on one half only: run stepAIC() on wn[half, ], take the terms with p < 0.01, refit a model with just those terms on the other half (wn[-half, ]), and report their p-values there. Which survive? What would you have reported from (a) alone, and what does this tell you about p-values from a selected model?

Hint

summary(s)$coefficients gives the p-values; reformulate(terms, "Flowering") builds the confirmation model. The noise columns are rnorm() draws: they have no relationship with flowering time by construction. Slides Pure noise, ‘significant’ results and Fix: choose on one half, test on the other; code page Part E, sections 35 and 36.

--- title: "Problem Set 5: Model Selection" subtitle: "PLS 206 — Applied Multivariate Modeling" author: "Grey Monroe" date: "Fall 2026" format: html: theme: flatly toc: true toc-depth: 2 toc-title: "Parts" highlight-style: github code-fold: false execute: echo: true eval: false --- ::: {.callout-note} **Due:** Monday, November 2, 1:10 PM (start of class) — submitted through Canvas. Questions 1–4 use Monday's lecture (AIC by hand, Akaike weights, AICc and BIC, testing a small set of hypotheses with `model.sel()`). Questions 5–10 use Wednesday's (collinearity and VIF, `dredge()`, importance and model averaging, stepwise selection and its traps). The [week's code page](Week5_Model_Selection.html) has working code for every step; the hints point to its parts and sections. ::: --- ## Submission Instructions Upload **two files** to Canvas: 1. **Written answers** — a `.pdf` document with your written responses, interpretations, and figures. 2. **R script** — a `.R` file with all the code you used to produce your answers. List the first and last names of any collaborators in both files. **Naming convention:** `HW5emailID.R` (code) and `HW5emailID.pdf` (written answers), where `emailID` is the part of your UC Davis email before the `@`. > **Example:** Grey Monroe would submit `HW5gmonroe.R` and `HW5gmonroe.pdf`. Feel free to use Google, StackOverflow, and other resources if you are stuck with code, and post questions on Slack if anything is unclear. **Every model you plot:** show the raw data, the fitted curve, and its 95% confidence band, with the response on the y-axis. --- ## The data: flowering time in wild wheat (again) You met these data in Problem Set 2. `Flowering` is flowering time (days from germination to flowering) for 54 genotypes of wild einkorn wheat (*Triticum monococcum*) collected across the Middle East and Europe, measured in a common garden: seeds from every collection site were grown together in a greenhouse, so differences among genotypes are genetic, not environmental. For each genotype's collection site we have `Latitude`, `Longitude` and the 19 bioclimatic variables of its home environment from [WorldClim](https://www.worldclim.org/data/bioclim.html) (`bio1` to `bio19`). In Problem Set 2 you tested one climate variable at a time. This week the question is which explanation is best supported, and what happens when we let the computer choose. ```{r} library(ggplot2) library(MuMIn) # model.sel(), dredge(), sw(), model.avg() library(car) # vif() library(MASS) # stepAIC() wheat <- read.csv("https://greymonroe.github.io/PLS_206/data/wild_wheat.csv") ``` --- ## Part 1: Three Hypotheses Three ideas for what sets flowering time in wild wheat, each with one variable to stand for it: - **Temperature:** warmer sites select for faster development (`bio1`, annual mean temperature). - **Summer drought:** plants from sites with dry summers flower early to escape the drought (`bio18`, precipitation of the warmest quarter). - **Day length:** photoperiod changes with latitude, and flowering is often cued by day length (`Latitude`). Your candidate set is six models: the null (`Flowering ~ 1`), each idea alone, temperature + drought (`bio1 + bio18`) and temperature + day length (`bio1 + Latitude`). **Q1.** Fit the temperature model, `lm(Flowering ~ bio1, data = wheat)`. Using `logLik()`, calculate its AIC by hand with AIC = −2 log L + 2k, and check your answer against `AIC()`. What is k for this model, and what does it count? ::: {.callout-tip collapse="true"} ## Hint `as.numeric(logLik(m))` gives log L, and `attr(logLik(m), "df")` gives the k that R uses. Count the parameters the model estimates: there is one more than the coefficients in `coef(m)`. Code page Part B, sections 6 and 7; slides *logLik() in R* and *AIC() does the same sum*. ::: **Q2.** Fit all six models and compare them with `model.sel()`. Which hypothesis is best supported, and what is its Akaike weight? The two-predictor models are both about 2 AICc units behind the temperature model: compare their log-likelihoods with the temperature model's and explain where the gap of about 2 comes from. Do drought or day length add anything once temperature is in the model? ::: {.callout-tip collapse="true"} ## Hint `model.sel(m0, m_temp, m_dry, m_day, m_temp_dry, m_temp_day)`. The `logLik` column is the fit, the `df` column is k. Slides *Only differences matter: ΔAIC*, *Akaike weights: share of the support*, *model.sel(): the six models, ranked* and *The evidence, as weights*; code page Part B, section 10, and Part C, sections 13 and 14. ::: **Q3.** Rank the same six models by AIC and by BIC as well as by AICc (`model.sel(..., rank = AIC)`, `rank = BIC`). Does the ranking change? How does the weight of the temperature model change from AIC to AICc to BIC, and why? With n = 54, which of AIC and AICc should you report? ::: {.callout-tip collapse="true"} ## Hint AICc = AIC + 2k(k + 1) / (n − k − 1); BIC = −2 log L + k log(n). Compare each one's penalty for one extra parameter: `log(54)`. Slides *AICc: a correction for small samples*, *BIC: a stiffer price* and *AIC or BIC?*; code page Part B, sections 11 and 12, and Part C, section 13 (`rank = BIC`). ::: **Q4.** Three comparisons. For each, say whether you can use `anova()`, AIC, both, or neither, and give the result where you can: - a) temperature alone vs temperature + drought; - b) temperature alone vs drought alone; - c) a labmate fits `lm(log(Flowering) ~ bio1, data = wheat)`, finds an AIC of about −146 against 353 for your temperature model, and says the log model is far better. ::: {.callout-tip collapse="true"} ## Hint `anova(m_temp, m_temp_dry)`. `anova()` needs one model to be a special case of the other. For c), ask what the likelihood is a likelihood *of*. Slides *Nested models: anova() and AIC agree*, *Not nested? Only AIC can compare them* and *Same response: log(y) and y cannot be compared*; code page Part C, sections 15, 16 and 18. ::: --- ## Part 2: Exploring All 19 Climate Variables **Q5.** Fit the model with all 19 climate variables. `bio7` gets `NA` (you found why in Problem Set 2: `bio7 = bio5 − bio6`), so drop it. Then prune the other 18 with `vif()`: drop the variable with the highest VIF, refit, and repeat until every VIF is below 5. Which variables are left? Was `bio1`, the best single predictor of flowering time, kept? Why did the pruning throw it out? ::: {.callout-tip collapse="true"} ## Hint A `repeat` loop saves you 12 rounds of copy-paste: ```{r} vars <- setdiff(paste0("bio", 1:19), "bio7") repeat { v <- vif(lm(reformulate(vars, "Flowering"), data = wheat)) if (max(v) < 5) break vars <- setdiff(vars, names(which.max(v))) } ``` For the last part, look at `cor(wheat$bio1, wheat[, c("bio5", "bio6", "bio10", "bio11")])`. Slide *Drop the worst VIF, refit, repeat*; code page Part D, section 25. ::: **Q6.** Run `dredge()` on the model with the variables left after Q5. What is the top model by AICc, and which models are within 2 AICc units of it? Rerun with `rank = "BIC"`: is the top model the same, and how does the runner-up differ? Plot `Flowering` against the top model's predictor with the smallest p-value, with the top model's fitted line and 95% band (hold the other predictors at their means). ::: {.callout-tip collapse="true"} ## Hint `options(na.action = "na.fail")` first, or `dredge()` refuses to run. `d <- dredge(global)`; `subset(d, delta < 2)`; `dredge(global, rank = "BIC")`; `top <- get.models(d, 1)[[1]]`. For the band, build `newdata` with that predictor on a grid and the others at `mean()`, then `predict(top, newdata, interval = "confidence")`. Slides *dredge(): every subset of the 7*, *The top 10 of 128* and *AICc vs BIC with 6,299 records*; code page Part D, sections 27 to 29, and Part B, section 6 for the band. ::: **Q7.** Report each variable's importance from `sw()` on the full dredge table, and average the models within 2 AICc units (`model.avg(d, subset = delta < 2)`). Every variable appears in exactly half of the 64 models; what importance would a variable that explains nothing get, and what does that say about the importance of `bio12`? In the model average, why is the `bio12` slope smaller in the full average than in the conditional average? ::: {.callout-tip collapse="true"} ## Hint `sw(d)`; `summary(model.avg(d, subset = delta < 2))` prints both averages. The full average counts a slope of zero for every model that leaves the variable out. Slides *Importance: summed weights* and *Model averaging: estimates across the top models*; code page Part D, sections 30 and 31. ::: **Q8.** Your Part 1 temperature model and the top dredge model from Q6 were fit to the same 54 genotypes with the same response, so their AICc values can be compared. Which is better supported, and by how much? Is the top dredge model a finding or a hypothesis, and which of the three goals (exploration, inference, prediction) does the dredge in Q6 serve? Two or three sentences. ::: {.callout-tip collapse="true"} ## Hint `AICc(m_temp)` and `AICc(top)` (both from `MuMIn`). Think about which variables `dredge()` was allowed to try. Slides *A tool for each goal* and *What we found: hypotheses, not results*; code page Part A, section 5, and Part D, section 32. ::: --- ## Part 3: The Traps **Q9.** Stepwise selection follows one path. (a) Run `stepAIC()` backward from the 18-variable model (all bio variables except `bio7`), and forward from the null model with the 18-variable model as the upper scope. Do they agree? Compare their AIC values. (b) Now take the six variables from Q5, resample the 54 genotypes with replacement 200 times (`set.seed(206)` first), rerun backward `stepAIC()` on each resample, and tabulate how often each variable is kept. Which variables are stable, and how often does the most common model come up? ::: {.callout-tip collapse="true"} ## Hint (a) `stepAIC(full18, trace = 0)` goes backward; `stepAIC(lm(Flowering ~ 1, data = wheat), scope = list(lower = ~ 1, upper = formula(full18)), direction = "forward", trace = 0)` goes forward. Compare with `AIC()`, not the AIC column `stepAIC()` prints (it drops a constant). (b) In a loop: `wb <- wheat[sample(nrow(wheat), replace = TRUE), ]`, refit and run `stepAIC()`, then `attr(terms(s), "term.labels")` gives the kept variables. Slides *Stepwise: one path through the maze* and *Resample the farms, rerun stepAIC*; code page Part E, sections 33 and 34. ::: **Q10.** Selection finds patterns in noise. Add 10 columns of pure noise to the six Q5 variables and draw a random half of the genotypes, running this block exactly as written: ```{r} set.seed(19) wn <- wheat[, c("Flowering", "bio3", "bio4", "bio6", "bio8", "bio12", "bio18")] for (j in 1:10) wn[[paste0("noise", j)]] <- rnorm(nrow(wn)) half <- sample(nrow(wn), 27) ``` (a) Run backward `stepAIC()` on `lm(Flowering ~ ., data = wn)` using all 54 rows. Is any noise column kept with p < 0.01? (b) Now select on one half only: run `stepAIC()` on `wn[half, ]`, take the terms with p < 0.01, refit a model with just those terms on the other half (`wn[-half, ]`), and report their p-values there. Which survive? What would you have reported from (a) alone, and what does this tell you about p-values from a selected model? ::: {.callout-tip collapse="true"} ## Hint `summary(s)$coefficients` gives the p-values; `reformulate(terms, "Flowering")` builds the confirmation model. The noise columns are `rnorm()` draws: they have no relationship with flowering time by construction. Slides *Pure noise, 'significant' results* and *Fix: choose on one half, test on the other*; code page Part E, sections 35 and 36. :::